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ABSTRACT 

The Fermi satellite has recently detected gamma ray emission from the central regions of our Galaxy. 

This may be evidence for dark matter particles, a major component of the standard cosmological 
model, annihilating to produce high-energy photons. We show that the observed signal may instead 
be generated by millisecond pulsars that formed in dense star clusters in the Galactic halo. Most 
of these clusters were ultimately disrupted by evaporation and gravitational tides, contributing to a 
spherical bulge of stars and stellar remnants. The gamma ray amplitude, angular distribution, and 
spectral signatures of this source may be predicted without free parameters, and are in remarkable 
agreement with the observations. These gamma rays are from fossil remains of dispersed clusters, 
telling the history of the Galactic bulge. 

Subject headings: 


1. INTRODUCTION 

While there are strong indications for the existence 
of cold dark matter from its gr avitational effects (e.g. 
iPlanck Collaboration ct ahll20l4 ) , there has not yet been 
any conclusive direct or indirect detection of the corre¬ 
sponding dark matter particles. One promising avenue to 
look for these particles is through annihilation in which 
two dark matter particles (a particle and its antiparticle) 
convert into high energy photons that we can observe. 
The dark matter annihilation signal is expected to be 
strongest where the density of dark matter is highest, 
i.e., in the centers of galaxies. 

Detailed analyses of the Fermi satellite’s map of the 
gamma-ray sky have revealed an excess around the 
Galactic center peaking at energies of GeV (e.g . 
Hooper fc G oodenoughl 2011 : IGordon fc Maciasl 120131 : 
Davlnn e: al . 201 II: . This excess appears to be roughly 
spherical and extends at least ~10-20° (1.5-3 kpc) 
from Sgr A*, the Galaxy’s central supermassive black 
hole. Remarkably, this signal can be interpreted as 
phot ons from annihilating ~30 GeV dark matter part i- 
cles (|Hooper fc Goodenoughl [201li Davlan et al.lf2014 i). 
In order to confirm this extraordinary interpretation, 
one must carefully rule out all other astrophysical 
sources. Possible alternatives include millisecond pul¬ 
sars (MSPs), rapidly spinning neutron stars that are ob¬ 
served in other regions of the Galaxy with very simi- 
lar gamma ray spectra to that of the observed excess 


(IGordon fc Macfasll2013l: Abazaiian & Kanlinghatl 

2012; 

Yuan & Zhand 12014 Abazaiian 

20111: Miraball 

2013 

Yuan & Iokai:2015l Petrovic et al.l 

20151): highly magne- 


tized young pulsars created in the innermost nuclear star 
clust er (lO’Learv et al.ll2015h ; injection of cosmic-ray pro- 
tons dGarlson fc Profumoll20l4 h or cosmic ray outbursts 
(iPetrovic et al.ll20l4 h However, it remains to be shown 
that any of these sources is sufficiently abundant and 
spatially extended to explain the gamma-ray excess. 

Energetic photons have also been observed from within 
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the central few pc around Sgr A* itself, extending from 
soft X - rays to ~100 Te V gam m a rays (jB agano ff et, al.l 
20011 lAharonian et all 12004 iBelanger et al.l 120061 
Perez et al. ns The origin of this emission is sub¬ 
ject to debate; see Ivan Eldikl (12015T) for a review. The re¬ 
gion near the event horizon of Sgr A* is likel y responsibl e 
for bright outbursts in soft X-rays (|Baganoff et al.ll200lt l, 
but this scenario struggles to explain the steady emission 
at much higher energies. Alternative explanations for the 
GeV and TeV flux inc lude the supernova remnant Sgr A 
East dCrocker et al.l [2005h . though this is strongly dis¬ 
favored based on its observed offset from the very hig h 
energy emission centered on Sgr A* (lAcero et al.ll2010 i. 
Secondary emission from particles accelerated by Sgr A* 
is another candidate, either in a steady state or from 
a past burst of accret i on (e . g. lAtovan fe Dermei 12004 
lAharonian fe Neronovl 120051 : iChernvakova et al.1 201 llj . 
Most of these scenarios cannot account for both the GeV 
and TeV emission. In contrast, a population of ^1000 
MSPs in the inner few pc coul d account for the emissio n 
from GeV through 100 TeV (|Bednarek fc Sobczakll2013l ). 
None of these scenarios seek to explain the GeV excess 
extending several kpc from Sgr A*. 


The pulsar population in the Galactic c enter has l ong 
been sought to test th e theory of gravity dPfahl fc Loebl 
12004 iLiu et ahl 120121 and references therein) and the 
existence of inte rmediate mass bla ck holes and grav¬ 
itational waves (iKocsis et al.112012T) . Extended multi¬ 
wavelength observations were conducted which should 
have detected a significant fraction of the most com¬ 
mon second-period pulsars, but only four were seen. 
This missing pulsar problem indicates that the for¬ 
mation and/or retention of ordinary pulsars may 
be i nefficient in this regi on (jDexter fe O’Learvl 12014 
iMacauart fe Ka nckar 2 015D . However, these searches did 
not significantly constrain the number of MSPs, espe¬ 
cially at the relatively large galactocentric distances of 
0.1-1 kpc where the gamma ray excess is observed. 

In this paper, we argue that the MSPs needed to 
produce the gamma ray excess were not made under 
the present conditions of the Galactic bulge, but were 
produced in dense globular clusters that have since dis- 
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solved. The population of globular clusters constitutes 
a key component in the theory of galaxy evolution, the 


dTremaine et al. 1974 Arca-Sedda & Canuzzo-Dolcettal 

2014 iCapuzzo-Dolcettal 19931: 

Lotz et al.l 20011; 

Bekki et al.l 2004 iCapuzzo-Dolcetta & IV 

occhi 

20081 Camizzo-Dolcetta & Mastrobuono-Battisti 

20Q9t 

Aearwal & Milosavlievic 2011b Hartmann et al. 

2ont 

Leieh et al.l 20121: Antonini et al. 


2012b Antonini! 

2014 20(3 den Brok et al. 2014 

Kruiissen 2014; 

Perets & Mastrobuono-Battistil 2014 Antonini et al.l 


have survived throughout the evolution of the Galaxy, 
and may be a small fraction of the initial cluster 
population. 

We organize the paper as follows. In Section [21 we dis¬ 
cuss the relationship of MSPs to globular clusters and the 
evolution of the Galaxy’s population of globular clusters. 
Section [3] discusses the scaling of the gamma ray lumi¬ 
nosity to the predicted population of disrupted clusters, 
while Section [4] presents the predictions of this model for 
the Fermi excess. Sections[5]and[6]discuss two objections 
to a MSP explanation of the excess, the high end of the 
luminosity function and the average spectrum. We dis¬ 
cuss prospects for radio detections of our predicted MSPs 
in Section [7] and conclude with Section [8] 


2. MILLISECOND PULSARS FROM DISRUPTED 
GLOBULAR CLUSTERS 

MSPs are thought to be “recycled” pulsars, spun up 
by t he accretion of material from a cl ose bi nary compan¬ 
ion dBhattacharva fe van den llouvellll99ll) . These close 
binaries are formed and driven to smaller separations in 
dense stellar environments where the rate of stellar dy¬ 
namical encounters is high. Once these interactions have 
sufficiently decreased the binary separation, the neu¬ 
tron star’s companion transfers material and angular mo¬ 
mentum, and reduces the neutron star’s magnetic field. 
This phase lasts ~10'-10 9 years and is visible, for low- 
mass companions, as a low mass X-ray binary (LMXB) 
(jlvanova et all [2008h . A long-lived MSP remains after 
the mass transfer stops. While the strong magnetic brak¬ 
ing in ordinary pulsars leads to a rapid sp i ndown, MSPs 
persi st for ^10 10 years dBhattacharva fe van den Heuvell 
EUD; these pulsars can long outlive their birth clusters. 

The highest abundances of MSPs are found in glob¬ 
ular clusters, the old dense stellar islands orbiting 
in the galactic halo. There are several indications 
that a fraction of the stellar mass of galactic bulges 
may have been formed by dissolving globular clusters 
dTremaine et all 1 19751 lArca-Sedda fe Capuzzo-Dolcettal 

12014 iGnedin et al.l120141 1 The distribution of old glob- 
ular clusters within galaxies and the geometry of galac¬ 
tic bulges are both approximately spherical. The num¬ 
ber density of globular clusters increases inwards on kpc 
scales, but shows a relative decrease within the galactic 
bulge. A tight correlation is observed between the mass 
of the galactic bulg e and the number of globular clusters 
dHarris et al.ll2014f ). 

Massive globular clusters spiral in towards the Galactic 
center due to dynamical friction. In the central kpc, the 
tidal gravitational field of the Galaxy may exceed the 
attractive field of the cluster stars, stripping the clus¬ 
ter from its outskirts and eventually down to its core. 


The cluster then spills its entire contents, including the 
MSPs in its core, into a spherical shell about the Galac¬ 
tic center. Since MSPs are long-lived, they remain bright 
in gamma rays after the cluster is disrupted. The high 
dynamical encounter rates needed to form new LMXBs 
and new MSPs are, however, strongly suppressed af¬ 
ter the disruption of the globular cluster. Therefore, 
the MSP population will be frozen at the time and or¬ 
bit of the cluster’s disruption while LMXBs (precursors 
to MSPs) will burn out within ~10 8 years. As a re¬ 
sult, the ratio of LMXBs to MSPs from disrupted glob¬ 
ular clusters becomes much lower in the galactic bulge 
than in the surviving globular clusters we observe to¬ 
day. Indeed, LMXBs are observed to be rar e in the bulge 
(iRevnivtsev et al.ll2008l iGholis et al.ll2015[ ). 

We model the distribut ion of globular c luste rs and the 
Galactic bulge following IGnedin et al.l d2014l) , who ac¬ 
count for mass loss from passive stellar evolution, clus¬ 
ter evaporation, infall due to dynamical friction, and 
tidal disruption; we adopt all of their fiducial parame¬ 
ters. This simple model was set up to reproduce the 
radial and mass distribution of extant globular clusters 
in the halo with no eye towards reproducing the Fermi 
excess. The distribution of mass from dissolved globu¬ 
lar clusters is shown in Figure 3 of lGnedin et al.l • 201 1 :: 
we use this result directly. It is a cored radial profile 
with p(r) ~ r ~ * 2 - 2 * and an enclosed mass ~10 8 M 0 at 1 
kpc. We assume that most gamma-ray sources in globu¬ 
lar clusters formed sufficiently quickly that the luminos¬ 
ity per unit mass in a population of disrupted clusters 
may be approximated by the observed value in surviving 
globular clusters. Their radial distribution is set by their 
orbits at their points of disrup tion i n the mo del. W e do 
not tune the fiducial model of lGnedin et al.l (120141 : our 
approach has no free parameters. 

The Appe ndix gives a brief overview of the 
IGnedin et al.1 (1 2014 11 model and calculations, with equa¬ 
tions giving the characteristic disruption timescales. We 

refer the reader to that paper for a more thorough dis¬ 
cussion. 

3. SCALING THE GAMMA RAY LUMINOSITY 

We compute the gamma ray luminosity p er unit stel¬ 
lar ma ss for the globular clusters studied by lAbdo et al.1 
(f20T(l and l isted in Table [T] Of these eleven clusters, 
lAbdo et~ahl (120101) reported gamma ray detections for 
eig ht. We use the mo re recent 2 GeV fluxes measured 
by iCholis et al.l (I2015D : this data set includes fluxes for 
two of the three globular clusters that were previously 
undetected . We ad opt the cluster distances compiled by 
lAbdo - et al.l (120 101 and convert the absolute V magni- 
tudes of lHarris! 19961 (2010 edition) to mass by assum¬ 
ing a mass-to-light ratio of 3 appropriate fo r 

an old, slightly metal-poor population (lMarastonll2005D . 
The mass of Terzan 5 is uncert ain due to its very large 
extinction (lLanzoni et al.l [2010f) . with those authors fa¬ 
voring a value 2 x 10 6 M 0 compared to the 3 x 1 0 5 
implied by the absolute V magnitude given by lHarrisI 
( 1996 ). but this has little effect on our results. 

Our approach gives a 2 GeV flux density at a distance 
of 8.3 kpc of 2 x 10 -15 GeV cm -2 s -1 Mq 1 . The variance 
on this value as estimated from bootstrap resampling is 
30%, with additional uncertainties from the cluster hum- 
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Projected Distance (kpc) 


TABLE 1 

Properties of Globular Clusters 


0.001 


0.01 


0.1 


Name 


Dist A io a 

(kpc) 


Dist H io b 

(kpc) 


M v h -F 2 GeV C 
(mag) 


47 Tuc 
uj Cen 
M 62 

NGC 6388 
Terzan 5 
NGC 6440 
M 28 

NGC 6652 
NGC 6541 
NGC 6752 
M 15 


4.0 ±0.4 

4.8 ±0.3 

6.6 ±0.5 

11.6 ±2.0 

5.5 ±0.9 

8.5 ±0.4 
5.1 ±0.5 
9.0 ±0.9 

6.9 ±0.7 
4.4 ±0.1 
10.3 ±0.4 


4.5 
5.2 
6.8 

9.9 

6.9 

8.5 

5.5 
10.0 

7.5 
4.0 
10.4 


-9.42 

-10.26 

-9.18 

-9.41 

-7.42* 

-8.75 

-8.16 

- 6.66 

-8.52 

-7.73 

-9.19 


5.6 

2.8 

3.8 

3.4 
12.6 

2.9 
3.8 

1.4 
0.9 
0.5 


a From 17\bdo ct all (120100 . references therein 
b From lHarrislll996l (2010 edition) 
c E 2 dN/dE , unit s are 10~ 9 GeVcm -2 s — h 
ICholis et al.l d2014D . 

* Uncertain due to very high extinction 


From 


nosities and distances, and possible systematic variations 
of cluster properties with galactocentric radius. Increas¬ 
ing the mass of Terzan 5 to 2 x 10 6 would d ecrea se this 
value by ~10%, while adopting the lHarrisI [19961 (2010 
edition) distances would increase it by a similar factor. 
There are also uncertainties about how representative 
these extant clusters are of the initial population and 
indeed about the initial population itself; we therefore 
(somewhat arbitrarily) adopt a factor of 2 (0.3 dex) as 
the uncertainty in our gamma ray luminosity scaling. 

We neglect systematic variations in the number of 
MSPs per unit globular cluster mass as a function of the 
cluster properties and evolutionary stage. The formation 



10 4 oo 


- 10 3 


Fig. 1.— Integrated flux (within an angle 4/ of the Galacti c 
center) of the Ferm i exc ess at 2 GeV by Gordon & Macias ( 20131) . 
ID avian et al.l (120141) and lAbazaiian et ah 1120145), compared to the 
prediction (solid blue curve) from disrupted globular clusters, as¬ 
suming the same gamma ray luminosity per unit mass as for intact 
clusters. The yellow hatching shows a factor of two uncertainty; 
the right axis shows the approximate number of enclosed MSPs. 
T he gamma ray signal from scaling the disrupted globular clusters 
of I Gnedin et al.l (120141) correctly predicts all meas urements, includ¬ 
ing an unresolved source around Sgr A* seen by [Abazaiian ct al. 
(|2014l) . with no free parameters. The black open stars include Sgr 
A* and the filled pentagons exclude it; we have also shown the 
disrupted globular cluster prediction with the inner 0?1 masked 
for comparison (blue dotted-dashed curve). We interpret this un- 
resolved sour c e as emiss ion from MSPs in a nuc lear star cluster. 
IDavlan et al.l (120141) and IGordon fe Maci'asl (120131') include this un¬ 
resolved source in their diffuse fits. 


0.5 


Projected Distance (kpc) 

1 1.5 2 2.5 


rate and the total number of MSPs are observed to cor¬ 
relate with the rate of encounters, not simply the cluster 
mass dHui et al.ll2010t iBahramian et al1l2013D , while the 
encounter rates are strongly affected by core collapse and 
the pri mord ia l binary fracti on which are not well under¬ 
stood (iBinncv fc Tremaine! I2008D . Indeed, the clusters 
listed in Table [T| have only a weak correlation between 
stellar mass and gamma ray luminosity. Future studies 
are needed to examine more detailed models of the MSPs 
within a population of globular clusters. 

4. THE PREDICTED FERMI EXCESS 

With an independent model of the p opulation of dis¬ 
rupted globular clusters (IGnedin et al.ll2014l ) and a scal¬ 
ing of total globular cluster mass to gamma ray lumi¬ 
nosity (Section [3j), we may compute the predicted Fermi 
GeV signal. Our results for the integrated flux within a 
circular region around the Galactic center, and the ap¬ 
proximate number of enclosed MSPs, are shown in Fig¬ 
ure G] as a function of the angular distance to the center. 
Figure[2]shows the differential flux within circular annuli. 
Each figure shows the recent measurements of the Fermi 
excess to be in excellent agreement with our predictions. 



Fig. 2.— Differ ential measuremen t s of t h e Galactic Center ex - 
cess at 2 GeV by Hoope r"fc SlatveH (120131) . IDavlan et al.l (120141) . 
and [Ualore ct al. (2 0151) , compared to the prediction from scaling 
the IGnedin et al.l (120141) disrupted globular clusters to the same 
gamma ray luminosity per unit mass as for intact clusters. We 
have included a factor of two uncertainty (yellow hatching) on the 
globular cluster prediction (blue curve). The blue curve and yellow 
hatching are not fitted to the data; they include no free parameters. 


The integral flux (Figure [Q contains two particularly 
noteworthy results. First, the red line in Figure [T| is 
a discontinous set of dark matter profiles fit at differ¬ 
ent annuli: the best annihilation fit to the data is not 
a physically meaningful dark matter profile. At sepa- 
ration s >0?5, which accou nt for nearly all of the flux 
in the IDavlan et al.l (|2014D fits, disrupted globular clus¬ 


ters provide a physical motivation for the form of this 
signal. Second, lAbazaiian et al.l (|2014D fit for an unre¬ 
solved source at the Galactic center, Sgr A*, which can 
be explained as MSP emission in a nuclear star clus¬ 
ter. We have ad ded this central s ource to the diffuse 
emission fitted by lAbazaiian et al.l (I2014D to obtain the 
open black stars in Figure [T] We have placed the left- 
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most black star at 0? 1, roughly the angular resolution of 
the Fermi telescope. The actual extent of the nuclear 
star cluster is predicted to be just a few pc, or ~0?02. 
To facilitate a direct comparison, we have also removed 
emission from the inner 0? 1 of the disrupted cluster pre¬ 
diction, and indicated the resulting diffuse gamma ray 
signal with a dotted-dashed blue line. Our prediction 
matches the gamma-ray flux at all separations. 

The existence of strong gam ma ray emission f rom a 
nuclear star cluster reported by lAbazaiian et ah . 120141 ) 
furt her supports the disrupted globu lar cluster hypoth¬ 
esis. iFaucher-Giguere fc Loebl (j2011[ ) found that the to¬ 
tal encounter rate in the inner parsec of the Galaxy is 
similar to that of the globular cluster Terzan 5, so the 
formation rate of MSPs (and resulting gamma ray lu¬ 
minosity) may be similar. The 2 GeV flux of Terzan 5 
would be just 5 x 10 -9 GeV cm -2 s -1 at 8.3 kpc. Un¬ 
less the gravitational potential of Sgr A* retains a much 
higher number of neutron stars than in globular clusters, 
this scenario can only explain ^10% of the flux seen by 
lAbazaiian et al.l (12014D and shown in Figure [TJ MSPs 
deposited in the nuclear star clus ter by massive g lobu- 
lar clusters were also suggested bv lBednarek fc Sobczakl 
(I2013H as an explanation of the observed TeV photons 
from around Sgr A*. Those authors required a popu¬ 
lation of ~1000-3000 MSPs, fully consistent with our 
results both in the Galactic center and at larger separa¬ 
tions. Recent 20-40 keV X-ray observations by the NuS- 
TAR satellite can also be expl ained by a large p opulation 
of MSPs in the inner few pc (jPerez et al.ll2015l ). 

5. THE MAXIMUM LUMINOSITY OF MSPS 

One objection to a population of bulge MSPs as the 
source of the Fermi excess is the paucity of individu¬ 
ally identified high luminosity gamma ray pulsars de¬ 
tect ed_as_poiiAsources within ~10° of t he Galactic cen¬ 
ter dCholis et alJl2015h . iLee et all (120151 1 found evidence 
that the GeV excess is from unresolved point sources 
with a 1.9-12 GeV flux cutoff at ~1.5-2 x 10 -1 ° pho¬ 
tons cm -2 s -1 , for a 0.1 100 GeV lumino s ity of ~1.5- 
2.5 x 10 34 ergs -1 at 8.3 kpc. iGholis ct all (12014 1 found 
a very hard luminosity function for MSPs, with most of 
the luminosity contributed by objects above a few 10 34 
ergs -1 , which should have been det ected as Ferm i point 
sources. We reanalyze the data of iCholis et al.l (I2014H 
a nd reexamine the cu toff at high luminosities. 

IGholis et all (|2014H use a sample of 59 field MSPs to 
determine the luminosity function (listed in their Table 
IV). The distance s they adopt are taken fr om the ATNF 
pulsar database (iManchcstcr ct ahl 120051 ) . available at 
http://www.atnf.csiro.au/p eople/pulsar/psrcat/ , 
which uses the model Galaxy of lTavlor fc Gordcsl ( 1993 1). 
hereafter TC93, to convert observed dispersion measures 
into distances. T he Fermi Seco nd Pulsar Catalog 
(2PC) distances (lAbdo et al.l I2013D listed in the same 
table generally use the same dispersion measures to 
calculate distances, but rely on the newer and more 
accur ate NE2001 model of the Galaxy dCordes fc Laziol 
120021 ). The newer distances are systematically lower, 
particularly out of the Galactic plane. Unlike the 
older TC93 model, the NE2001 model supplies enough 
free electrons to account for the obs erved dispersion 
meas u res of nearly a ll Gal actic pulsars dCordes fc Laziol 
I2002D . ICholis et al.l (|2014l) only used the 44 MSPs with 


TABLE 2 

Distances to Field Millisecond Pulsars 


Name 

DM 

pc cm -3 

£^093 

(kpc) 

t?NE2001 

(kpc) 

Ref 3, 

J0307+7443 

6.4 

0.34 

0.6 

R12 

J0533+6759 

57.4 

6.66 

2.4 

R12 

J0605+3757 

21.0 

1.16 

0.7 

R12 

J1137+7528 

29.2 

19.53 

1.5 

* 

J1142+0119 

19.2 

2.04 

0.9 

R12 

J1301+0833 

13.2 

0.91 

0.7 

R12 

J1302—3258 

26.2 

1.86 

1.0 

R12 

J1311—3430 

37.8 

3.72 

1.4 

R13 

J1312+0051 

15.3 

1.15 

0.8 

R12 

J1543—5149 

50.9 

1.46 

2.4 

N14 

J1544+4936 

23.2 

2.30 

1.2 

R12 

J1630+3734 

14.1 

0.85 

0.9 

R12 

J1640+2224 

18.4 

1.15 

1.16 

L05 

J1732—5049 

56.8 

1.81 

1.3 

V09 

J1745+1017 

23.9 

1.36 

1.3 

R12 

J1811—2405 

60.6 

1.70 

1.8 

N14 

J1816+4510 

38.9 

4.20 

2.4 

R12 

J1843—1113 

60.0 

1.97 

1.7 

H04 

J2129—0429 

16.9 

1.03 

0.9 

R12 

J2256—1024 

14 

0.91 

0.65 

B13 


a Re ferences abbreviated as: B13 IBrcton ct al. 

2013) ; C14 UCholis^et^lJJ^OlJ); H04 I Hobbs_eW4J 
2004); L05 IILo hmer et al.| 120051) ; N14 l|Ng et al.l 

2014) : R12 IRav et al.ll2012l) : R13 IIRav et alj|2013) : 


V09 IIVerbicsrct^nr2009l). 


http://astro.phys.wvu.edu/GalacticMSPs/GalacticMSPs.txt 

see text tor details 


|6| > 10° to determine the luminosity function; the 
revised distances thus have a large effect on the results. 

Table [2] lists the NE2001 distances for al l pulsar s 
without 2PC distances listed in ICholis ct all (I2014D . 

In the case of J1843—1113, iHobbs et all (j2004l ) have 
provided a distance using the TC93 model Galaxy. We 
have taken their dispersion measure and converted it 
into the tabulated NE2001 distance estimate. One MSP, 
J1137+7528, does not appear in any database. The only 
reference to this pulsar that we were able t o find (apart 
f rom the Fermi Third Point Source catalog, lAccro et all 
[20T5l) was in a table of pulsar data available at 
http :// astro.phys. wvu .edu/GalacticMSPs/GalacticMSPs.txt 
This site lists a dispersion measure of 29.2 pc cm -13 
which, according to the NE2001 free electron model, 
implies a distance of 1.5 kpc. 

We also note that these new distances have signifi¬ 
cant uncertainties. The brightest MSP in our sample, 
for exam ple, is .106 14—3329, w hich was identified with 
Fermi by iRansom et all (|2011l ) . This is a generally un¬ 
remarkable MSP apart from the very high gamma ray 
efficiency (greater than unity) im plied by i ts dis p ersion 
measure distance of 1.9 kpc. As IRansom et all ( 1201 1) 
note, J0614 is nearly tangent to the Gum nebula where 
the NE2001 model has a very steep gradient in dispersion 
measure (and hence, in derived distance); those authors 
suggest that the true distance is likely to be closer by a 
factor of 2 or more. Decreasing its distance by such a 
factor would reduce its integrated gamma ray luminosity 
to ~10 34 ergs -1 and remove the need to invoke a very 
large gamma ray efficiency and/or beaming factor. 

Other pulsars, including those just below the cutoff, 
could also have distance errors. Random errors, how¬ 
ever, will tend to smear out a distribution; a deconvo¬ 
lution should sharpen the cutoff. Also, an anomalously 
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Fig. 3.— Cumulative number of MSPs in thelCho 
sample with |6| > 10° adopting the Taylor &; Cordei 


blue line) and an updated model ( Cordes £; LazicT 


isetal ]4m3) 
(Iggl) (TC93, 


20021 NE2001. 


red line) of the Galaxy’s dispersion measure. The latter leaves just 
two MSPs sufficiently luminous to detect as point sources near the 
Galactic center. One, JQ614—33 29 is likely to be at l east a factor 
of 4 less luminous than shown (I Ransom et al.M201U) ; the other, 
J0218+4232, is just marginally (^ lcr) more luminou s than the 
Galactic center detection threshold JAbdo et al.l 120131) . The ver¬ 
tical shading shows the luminosity cutoff needed to repr oduce the 
gamma ray excess w i th an unresolved source population (|Lee et al.l 
120151 ') : Bartels et al. (2019) favor a somewhat higher cutoff. 


large distance (by a factor of 2 , say) requires simply an 
unmodeled clump of ionized gas along the line-of-sight. 
An anomalously small distance requires a void or bubble; 
the required void for the same fractional distance error 
becomes larger with increasing distance. A factor of 2 er¬ 
ror in a 1 kpc distance would require a 1 kpc completely 
empty void (to go along with the 1 kpc of properly mod¬ 
eled free electron density). 

Figure [3] s hows t he cutoff at high luminosities in the 
iGholis et al.l (12015H MSP sample using both the TC93 
and the NE2001 distances. For reference, we have 
als o shown t he ap proximate luminosity cutoff favored 
bv lLce et alJ (12015H of ~1.5-2 photons cm -2 s -1 between 
1.9 and 12 GeV, transformed into 0.1-100 GeV flux as¬ 
suming the well-measured spectrum of J0614—3329 and 
placed at a distance of 8.3 kpc. There are only two 
MSPs above this cutoff. One, J0614 —3329 itself, was dis¬ 
cussed above; iRansom et al.l (|201lD argue that its true 
distance is likely to be at least a factor of ~2 smaller 
than that implied by its dispersion measur e. The other 
is J0218+4232, for which lAbdo et all (120131 1 quote a fac¬ 
tor of 2 error in the luminosity as a result of dispersion 
measure uncertainties in the NE2001 model, making it 
more luminous than 2 x 10 34 ergs -1 by just ~1 <t. The 
updated puls ar distances sup port the luminosity cut ¬ 
off needed by iLee et all (120151) and iBartels et al.l (|2015H 
to accou nt for the GeV ex cess with unresolved point 
sources. IBartels et alJ (120151) fav or a cutoff at higher lu- 
minosities than lLee et al.l (|2015t l (though over a different 
bandpass). 

6 . THE AVERAGE SPECTRUM OF MSPS 

We now turn to the spectrum of an unresolved distribu¬ 
tion of Fermi MSPs, esti matin g the int egrate d light using 
the field MSPs listed in ICholis et alJ (I2014H . The spec¬ 
trum of MSPs has been suggested to differ modestly from 


that of the observed GeV excess (ICholis et al.ll2015lf , ar¬ 
guing against their ability to produce the gamma rays 
around the Galactic center. We use the field MSPs under 
the assumption that many of them have similar origins 
to the one we suggest for the central bulge population, in 
the cores of stellar clusters disrupted long ago. We also 
account for the fact that Fermi s sensitivit y is a strong 
function of photon energy (lAcero et al.ll2015f) . Following 
ICholis et al.l (12014D . we do not apply a cut in Galactic 
latitude b (as we and they did for the high luminosity 
cutoff, Section 0. Applying such a cut would make the 
spectra we show in t his section s lightly harder, and lessen 
the tension with the lDavlan et alJ(|2014f) spectrum of the 
GeV excess. 

The MSP sample of lCholis et al.l (12014H is not selected 
at a single frequency: for a fixed 2 GeV flux density 
with little emission from higher energy photons, very 
soft sources are easier to detect than harder sources. 
Only the soft sources supply enough photons at energies 
<1 GeV to contribute significantly to their detectabil¬ 
ity. We therefor e creat e a 2 GeV-selected sample from 
the ICholis et al.l (|2014f ) MSPs by including only those 
with a 1-3 GeV test statistic of 8 2 , corresponding to a 
signal-to-noise ratio of at least 8 in this band. Of the 59 
MSPs, 45 pass this cut. This sample of 45 MSPs should 
be complete independently of spectral shape. 

Under the (unlikely) assumption that our field MSPs 
form a flux-limited sample of a spatially uniform pop¬ 
ulation, we can derive a weighted average of the in¬ 
dividual MSP spectra that matches the spectrum ex¬ 
pected for a population at fixed distance. Uniformly 
weighting each MSP’s spectrum is equivalent to assum¬ 
ing a volume-limited survey. These two scenarios, flux- 
limited/spatially uniform and volume-limited, almost 
certainly bracket the truth. The rest of our calculations 
use an average or composite of these two limiting cases. 

Assuming a flux-limited survey, the total gamma-ray 
flux from Fermi -resolved MSPs in a luminosity range 
[L, L + SL] is 

„ 2 L dN dN 

F tot oc / drr——dL = r ma ^L—dL, (1) 
J 0 r z dL dL 

where dN/dL is the differential number density and 
r max oc \J L/F m i n is the maximum radius out to which 
the source could be detected. For a source population 
around the Galactic center, we wish to know the total 
luminosity per unit volume, which is simply given by 

L dN „ 

Fqc oc - 5 —— SL . ( 2 ) 

r GG dL 

We combine Equations 0 and (0 to obtain 

Fee oc F to tL~ 1/2 . (3) 

The integrated spectrum from an unresolved MSP pop¬ 
ulation near the Galactic center may therefore be esti¬ 
mated by scaling the observed flux density of each ob¬ 
ject by that object’s luminosity to the — 1/2 power, and 
simply adding all of the scaled flux densities together. 

Figure 0 shows the results, with all spectra normalized 
to their 1.7 GeV values. The blue dotted-dashed line 
shows the result for simply adding together all MSP spec¬ 
tra without selecting them by 1-3 GeV flux and without 
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scaling by luminosity. The red line is the average of the 
spectra with and without weighting by L -1 / 2 , i.e, as¬ 
suming volume-limited and flux-limited samples, respec¬ 
tively. The blue and orange hatching show the ler and 
2 cr uncertainties in the red spectrum as estimated from 
bootstrap resampling of the 45 MSPs. For this exer¬ 
cise, we have adop ted the fitted spectra in Table I of 
iCholis et al.l (|2014f) and have neglected measurement er¬ 
rors, fitting errors, and distance errors. 

The difference in Figure [4] between the scaled and un¬ 
sealed spectra results from a correlation between lumi¬ 
nosity and spectral index. Distance errors will tend to 
blur this correlation; the MSP spectrum of a population 
at a single distance is likely to be slightly harder than the 
red line in Figure^ Including this effect and adding mea¬ 
surement errors would not bring the MSP spectrum into 
perfect agreement with the Galactic center excess, but it 
could bring the ler discrepancy to as little as ~20-30% 
at 500 MeV. Selecting only those MSPs with |6| > 10° 
(38 of the 45 that pass our 1-3 GeV signal-to-noise cut) 
would also marginally improve the agreement with the 
spectrum of the GeV excess. 

The discrepancy between our estimated average MSP 
spectrum and the GeV excess is only significant at the 
lowest energies (<800 MeV) where Fermi 's sensitivity is 
rapidly falling. Unce rtainties in G alacti c diffuse emis¬ 
sion are largest here dCalore et al.l 120151) . As a result, 
there are spectrally correlated systematic errors in the 
spectrum of the GeV excess not shown in the black 
stars of Figure [4] Systematic errors can be quite large, 
and can also arise from the method of masking point 
sources and from the assumed morp hology of the excess , 
among other aspe cts of the fitting (jDavlan et al.l 120141 : 
iCalore et al.ll2015h . Figure 0] also shows the systematic 
err ors from vary i ng the diffuse backgrounds as estimated 
by I Galore et al.l ( 2015 1. These gray and gold hatched 
regions neglect statistical errors. 

7. PROSPECTS FOR RADIO DETECTIONS 

Our results show that a population of disrupted glob¬ 
ular clusters, which must exist to explain the current 
clusters, naturally predicts a field population of MSPs in 
the Galaxy’s inner few kpc. These MSPs satisfy the spa¬ 
tial, spectral, and luminosity requirements imposed by 
the Fermi observations. A large population of MSPs in 
a nuclear star cluster is another necessary consequence 
of a population of disrupted massive globular clusters. 
Such a population expla ins the 20-40 keV X-ray emis¬ 
sion seen by NuSTAR (jPerez et al.l 120151) and implies 
that many of the unidentified Chandra point sou rces may 
be MSPs dMuno et al.l 120041 : iPerez etahl 120151) . Astro- 
Id (jTakahashi et alJ 120101 ) will also be sensitive to high- 
energy X-rays, and could confirm the NuSTAR results. 
A population of ~1000 MSPs around Sgr A* can also 
explain the observed TeV emission by inverse Comp- 
ton scattering of the dens e interstellar radiation field 
(iBednarek fe SobczaklteORll) . 

Radio observations could individually detect our pre¬ 
dicted MSPs and confirm their identities. However, 
the bulk of the radio observations to date have fo¬ 
cused not on scales of tens to thousands of pc, where 
most of our predicted MSPs lie, but in the inner¬ 
most pc. This was motivated by theoretical estimates 
predicting ~100-1000 pulsars formed in situ within 
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Fig. 4. — The average spectrum of F ermi-detected field MSPs 
adopting the fitted spectral parameters of ICholis et al.l <120141) . The 
dotted-dashed blue line is the unweighted average spectrum. The 
red line has selected only those MSPs detectable based only on 
their 1-3 GeV flux (45 of 59 MSPs), and is the average of the spec- 
tra expected for a population at uniform distance assuming the 
ICholis et al.l (120140 to be volume-limited and flux-limited. These 
scenarios almost certainly bracket the truth. The blue and or¬ 
ange hatching show ler and 2cr sample variances as estimated us¬ 
ing bootstrap resampling. We have neglected errors in the MSP 
distances and in the spectral measurements; both would tend to 
alleviate the discrep ancy with the observed Galactic center excess 
JDavlan et al.l 12014) . The error bars on the IDavlan et al.l (120141) 
fits are only statistical; systematic errors (which are spectrally cor¬ 
related) are neglected. The gold and gray hatching show ler and 2cr 
sy stematic unce r tainti es (neglecting statistical errors) as estimated 
bv ICalore et al.l (|2015l ). 


0.02 p c of Sgr A* dPfahl fc Loebl 12004]) . More re¬ 
cently, iFaucher-Giguere fc Loebl (120111 ) noted that the 
encounter rate in the inner 1 pc of the central star clus¬ 
ter is comparable to that of the globular cluster Terzan 5 
(which has many MSPs), and estimated that up to ^1200 
MSPs may be present in this region due to the deeper 
gravitational potential well of Sgr A*. The disrupted 
globular cluster scenario instead predicts these MSPs to 
be found over a larger region: we predict ~1,000 MSPs 
within 3 pc of Sgr A*, and a further ~1,000 MSPs within 
300 pc (2°, see Figure [T]). 

MSP observations towards the Galactic center are ex¬ 
tremely challenging because of the large dispersion mea¬ 
sures. Radio pulses at a frequency v are broadened by 
an amou nt r = (1.3 ± 0.2 )(i//GHz)- 3 - 8±0 - 2 (with r in 
seconds, iSpitler et al.l 12014 ). implying that MSPs may 
not be observed below ~8 GHz. The radio intensity of 
pulsars scales steep ly with frequency (/ oc v~ ls to is~ l s , 
IKramer et al.lll998h . so high-frequency detections require 
extended integration times. 

While discovering and timing MSPs 0.001 pc from 
the central supermassive black hole would offer tanta¬ 
lizing measurements of general relativity and tests of 
alternative theories of gravity dWex fc Kopeikinl 1 1 9991 
IKramer et al.l 12004 iCordes et alJ 12004 iPfahl fc Loebl 
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12004 iLiu et all 120121 ). discovering MSPs further out 
within 10 pc would also be invaluable. Such MSPs 
could be used to measure the properties of the nu¬ 
clear star cluster, find intermediate mass black holes, 
and me asure the gravitat ional waves of the Galactic 
center (lKocsis et al.l l2012f ). If the nuclear star clus¬ 
ter indeed formed from disrupted globular clusters, we 
predict ~1000 MSPs within ~10 pc of the Galac¬ 
tic center; such a population of MS Ps can account 
for th e unresolved Fermi flux seen by lAbazaiian et all 
: 201 1Future high frequency radio surveys will have 
the frequency coverage a nd sensi tivit y needed to de¬ 
tect this MSP population dChennamangalam fe Lorimed 
l20l4lMacauart fc Kanekarll2015l) . Even larger radio sur- 

veys such as the square kilometer array (SKA) on 100 pc 
to 2 kpc scales are required to confirm or disprove the 
disrupted globular cluster origin of MSPs in the Galactic 
bulge. 

8 . CONCLUSIONS 

The Fermi Galactic center excess is in excellent agree¬ 
ment with independent predictions of the population of 
MSPs produced in disrupted globular clusters. This as- 
trophysical model appears to fit the observations as well 
as dark matter annihilation, but without any free pa¬ 
rameters. MSPs from disrupted clusters also provide an 
excellent match to the observed emission near Sgr A* 
from hard X-rays through very hard gamma rays. If the 
bulge indeed contains a large population of stars from 
long-dead clusters, such MSPs form a background that 
must necessarily be present in the Fermi data. 

The observed emis sion extends at least ~2 kpc from 
the Galactic center dHooper fc Slatverl 120131) . far from 
the nuclear star cluster around Sgr A* where dynami¬ 
cal formation of MSPs is plausible. LMXBs burn out 
after the disruption of globular clusters, reducing their 
relative numbers in the galactic bulge, consistent with 
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the lack of LMXB observations (c.f. iCholis ct al.ll2015 ). 
We conclude that the dominant MSP population is not 
likely to have formed under the current conditions in the 
bulge, but was deposited by dissolving globular clusters. 
If the Fermi excess is indeed the relic of a previous large 
population of globular clusters, it provides the first di¬ 
rect evidence for their existence, and strongly supports 
the theory for the globular cluster origin of the nuclear 
star cluster. Future radio observations may be directly 
sensitive to these MSPs and could offer decisive evidence 
of a broad distribution of MSPs deposited by globular 
clusters. 

While our results disfavor a dark matter interpretation 
of the GeV excess, they show that Fermi can offer a new 
probe of the formation history of the bulge, and of the 
evolution of the Galaxy’s globular cluster system. Our 
reevaluation of field MSP lum ino sities, combined wit h 
the results of iLee et~ahl (120151) and iBartels et al.l (12015D . 
suggest that we will soon begin to resolve the brightest 
of these fossils. 
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APPENDIX 


THE EVOLUTION OF A POPULATION OF GLOBULAR CLUSTERS 

We adopt the semianalytical method of IGncdin et al.l (I2014T) to calculate the evolution of a population of globular 
clusters in the Milky Way. The only result we use is the mass distribution of disr upted clusters , which we take directly 
from their Figure 3. We provide a brief summary of the basic assumptions of the lGnedin et al.l (120141 ) model here, and 
refer the reader to that paper for details. Using parameters that match the observed properties of young star clusters, 
the model recovers the present-day observed masses and radial distribution of globular clusters, and it predicts the 
total radial mass distribution of globular cluster mass deposited in the Galactic bulge. The initial mass distribution 
of globular clusters is set to dN/dM oc M~ 2 between 10 4 10 7 Mq, and the initial radial distribution of GCs is set to 
follow the mass distribution of Galactic field stars scaled down uniformly to 1.2% assuming that all globular clusters 
formed at z = 3. This distribution places equal mass in logarithmic bins, i.e., there is the same amount of stellar 
mass in 10 6 - 10 7 M 0 globular clusters as in 10 4 10 5 Mq clusters. All components (globular clusters, field stars and 
dark matter) are spherically distributed around the Galactic center; the field stars follow a Sersic profile with an 
index n s = 2.2, the enclosed mass is approximately M(R) ~ 10 5 Mq(R/100 pc) 2 ’ 4 if i? < 1 kpc; the dark matter 
follows an NFW profile pdm ex (i?/i? s ) -1 ( 1 + R/R s )~ 2 where R s = 20 kpc and a total mass within 12 R s is 10 12 Mg 
dNavarro et al.l[T~997T ) : and there is a supermassive black hole of 4 x 10 6 Mq at the center. 

Once the model is initialized, the clusters move on circular orbits in the instantaneous gravitational field of dark 
matter, stars, globular clusters, and the deposited mass fr om glob ular clusters, an d they spiral inwards due to dynamical 
friction. The inspiral time is proportional to M -1 (iBinnev fe Tremaind 120081) . implying that more massive clusters 
segregate inwards more quickly. The globular clusters evolve due to mass loss through stellar evolution, they slowly 
evaporate independently of their location relative to the galactic center, and they lose mass through tidal stripping by 
the galaxy. The orbital time, dynamical friction time, isolated cluster evaporation time, and tidal disruption time are 
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given, respectively, as 


torb = 0.06 Gyr- 


R / v c (R,M) 


kpc \100kms _1 


10 5 Mq 


-l 


t iso = 17 Gyr 
ttid = 67 


2 x 10 5 Mq 
to 


2 x 10 5 Mq 


torb{R , M) 


(Al) 

(A2) 

(A3) 

(A4) 


where R is the orbital radius, m is the cluster mass, M is the enclosed mass, a = 2/3 and f e = 0.5, v c (R,M) = 
( GM/R) 1 / 2 is the circular velocity. The clusters are followed individually and modeled using their average properties: 
half-mass radius, total mass, and average density, as 


dm to 

dt min(f t id,£iso ; twind) 

dR 2 _ R 2 
dt tdt 

Here t w ind is the timescale on which mass is lost due to stellar evolution and winds following dPrieto fe Gnediiill2008f) . 
Approximately 30% of mass is lost during the first 0.3 Gyr, and another 10% during the following 10 Gyr. The 
circular velocity and enclosed mass are updated as the radial mass distribution changes during the inward migration 

and evaporation of globular clusters. __ 

The average density of the cluster is assumed to vary with mass as (jGnedin et al.l [2014h : 


(A5) 

(A6) 


Mq 


Ph = 10' „ 

pc J 


rh = 


x min 1^10 , max 

/ 3to V /3 
\8np h (m) 


10 5 Mq 


(A7) 

(A8) 


A globular cluster is disrupted when the mean enclosed density in 
density of the cluster itself, 


c v c (R,M)) 2 
Ph 2 -kGR 2 


dark matter, 


gas, and stars exceeds the average 


(A9) 


Heavy, dense clusters can sink closer to the center before getting disrupted. The mass weighted mean lifetime of 
disrupted clusters is typically several Gyr. Thus, an average disrupted cluster may have had a similar number of MSPs 
per unit mass as similar mass clusters further out which survived disruption until the present (see Section O. 

In this model, the surviving globular clusters have an approximately lognormal mass distribution and a nonuniform 
radial distribution that is consistent with observations. The globular clusters that do not survive are typically disrupted 
before they reach the very center of the Galaxy, creating a characteristic cored density profile. The mass of the disrupted 
globular clusters exceeds the initial stellar mass in the the nuclear star cluster, the very central region of the galactic 
bulge, deliverin g a few 10 1 Mq wi thin ~10 pc of Sgr A*. An additional ~10 8 Mq is deposited interior to ~1 kpc 
(see Figure 3 in lGnedin et al.ll2014f) . The mass from disrupted clusters is deposited roughly spherically with a density 
decreasing with radius approximately as p ~ r~ 2 2 at 1 kpc. Here, the exponent depends on the assumed details of 
the initial cluster population, but is roughly constant (within a few tenths) between ~200 pc and a few kpc (~1° and 
20° projected at 8.3 kpc). 

While this toy model captures many of the essential features of globular cluster evolution within galaxies, it neglects 
several possibly important details. These include core collapse; binary and multibody interactions which may heat the 
cluster or eject stars; gas effects (accretion, inflow, star formation); resonant interactions, violent relaxation, radiative 
or thermal feedback from supernova explosions or an active galactic nucleus; the effects of galactic anisotropy (disk, 
bar, spiral arms); the effect of tidal shocks when crossing vertically through the galactic disk; the collision of globular 
clusters; the formation of new star clusters; the effects of galaxy mergers; and supermassive black hole binaries. 
Some clusters do display indications of interesting format i on histories that are hinted at by these complicating effects 
but are not captured in the toy mod el ()Bedin et al.ll2004 iMarino et al.l 120081 : iFerraro et al.ll2009l : IMarino et al.1120091 : 
lYong et al.l I20l4 iMilone et al.ll2()i~5ll . Some of these effects may have an influence on the radial mass-loss profile, 
but most will not affect the predicted spherical morphology of the tidal debris, as long as the initial distribution of 
globular clusters is roughly spherical and the rate of mass loss is slow over an orbital time. For example, a possible 
source of asphericity may result from tidal shocks generated by the cluster crossing the Galactic disk, which catalyses 
evaporation or core collapse. However, the characteristic timescale of tidal shocks ranges between 3 to 10 half-mass 
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relaxati on ti mes dGnedin et al.|[l999l) . which is between 0.1 and several Gyr, longer than the orbital time within 2 kpc 
(see Eq. [All) . The tidal debr is from g l obular clus ters overlapping with the gamma-ray excess may remain spherical in 
a wider class of models than lGncdin et ahl (|2014f) . 








